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I.  INTRODUCTION 


A,  OVERVIEW 

Since  the  start  of  launehing  vehieles  into  spaee,  there  has  been  an  ongoing  effort 
to  reduee  the  eost,  safety,  and  reliability  of  a  reusable  launch  vehicle  (RLV).  One  aspeet 
that  this  research  is  looking  to  correct  is  the  time  required  to  develop  optimal  launeh 
trajectories.  Optimal  launch  trajectories  are  essential  to  ensure  that  the  most  eost  effeetive 
launeh  trajeetory  is  flown.  In  this  research,  the  algorithm  that  will  be  used  is  DIDO. 
DIDO  is  a  MATLAB  optimal  eontrol  toolbox  that  was  named  after  Dido,  the  founder  and 
first  queen  of  Carthage.  She  is  famous  for  her  use  of  mathematies  in  solving  an  optimal 
eontrol  problem  (OCP)  before  ealeulus  was  even  invented.  DIDO  is  based  on 
pseudospeetral  optimal  eontrol  theory  that  is  designed  to  solve  an  OCP  in  the  same 
manner  as  using  equations  on  a  pieee  of  paper  [1].  The  diffieulties  in  solving  for  costates 
are  eliminated  by  the  eonveetor  mapping  prineiple  therefore  DIDO  produces  speetrally 
aceurate  solutions  [2].  With  this  tool,  a  more  eonvenient  method  to  determine  launch 
trajectories  ean  be  developed  to  help  reduce  the  time  spent  on  the  solution  of  a  launeh 
trajeetory. 

There  is  a  great  demand  for  satellite  based  equipment  and  it  only  keeps  getting 
larger.  The  military  is  heavily  reliant  on  launeh  vehieles  sinee  a  vast  majority  of  its  net- 
eentrie  warfare  is  based  on  satellite  eommunieations  [3].  A  great  deal  of  U.S.  national 
security  surveillance  is  done  via  satellite  for  the  ability  to  gather  the  most  real  time 
informational  available  [4].  The  global  positioning  system  (GPS)  is  not  only  vital  to  the 
military  but  to  the  eivilian  realm  as  well.  The  eivil  maritime  and  aviation  eommunities 
rely  heavily  on  GPS  for  aecurate  positioning  for  reliability  and  cost  savings.  Lastly, 
NASA  has  a  huge  demand  for  launeh  vehieles  as  they  are  responsible  for  resupplying  the 
International  Spaee  Station  (ISS)  and  sending  probes  for  deep  spaee  exploration,  as  well 
as  other  satellite  missions  sueh  as  the  James  Webb  Spaee  Teleseope  [5]. 
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B, 


DIFFICULTIES  WITH  THE  CURRENT  APPROACH 


The  current  cost  to  launch  a  pound  of  payload  into  an  Earth  orbit  is  around 
$10,000  [6].  In  order  to  reduce  the  total  cost  of  launch,  industry  is  constantly  looking  to 
reduce  the  mass  of  the  objects  being  sent  into  orbit.  The  other  aspect  is  to  reduce  the  cost 
to  launch  that  object  into  space.  This  is  the  primary  reason  why  companies  are  trying  to 
optimize  launch  vehicle  trajectories.  The  current  industry  standard  for  optimizing  launch 
trajectories  is  the  NASA  program  to  optimize  simulated  trajectories  (POST)  [7].  This  is 
an  immensely  complicated  program  that  takes  months  to  understand  how  to  operate. 
POST  takes  the  position  of  using  a  direct  shooting  method  to  calculate  state  variables  as  a 
function  of  time  [8].  Another  aspect  is  that  POST  requires  an  initial  guess  for  each 
independent  variable  that  would  otherwise  be  held  constant.  Developing  the  initial  guess 
can  be  very  time  consuming  [9].  That  complication  leads  to  how  long  it  takes  to  develop 
a  launch  trajectory  and  the  intense  man  power  required.  A  successful  launch  would 
require  being  able  to  predict  conditions  months  in  advance.  If  launch  conditions  are 
outside  of  those  that  were  predicted,  the  launch  may  have  to  be  terminated. 

C.  OBJECTIVE 

This  thesis  research  was  done  to  target  the  method  in  which  launch  trajectories  are 
developed.  The  goal  is  to  use  a  modem  algorithm,  DIDO,  to  reduce  the  time  that  is 
required  to  develop  launch  trajectories.  DIDO  removes  the  traditional  shooting  method  to 
solve  the  OCP  by  using  pseudospectral  optimal  control  theory  [1].  By  being  able  to 
develop  a  trajectory  closer  to  the  launch  date  allows  for  a  more  accurate  prediction  of 
conditions  to  develop  a  more  accurate  trajectory.  This  will  drastically  reduce  the 
manpower  costs  to  become  trained  on  the  software  and  develop  trajectories. 

Another  aspect  that  this  thesis  research  contributes  to  is  a  move  towards  more 
automation.  The  goal  being  that  the  algorithm  is  robust  enough  that  the  only  portions  that 
need  to  be  changed  are  the  starting  conditions  and  the  endpoint  conditions.  This  further 
increases  the  simplicity  of  the  method  to  solve  the  given  problem. 

This  thesis  research  specifically  addresses  the  goal  of  maximizing  the  first  stage 
final  velocity.  This  problem  was  chosen  in  an  effort  to  obtain  a  launch  vehicle  final 

2 


velocity  that  was  closer  to  the  final  orbital  velocity  in  a  more  expeditious  manner.  To 
account  for  real  world  uncertainties,  a  Gaussian  process  was  used  in  a  Monte  Carlo 
simulation  to  allow  the  worst  case  performance  to  be  identified.  This  knowledge  will  lead 
to  more  flexibility  in  the  launch  window  and  a  more  reliable  launch  trajectory. 

D,  THESIS  OUTLINE 

This  thesis  is  written  in  a  manner  that  shows  the  reader  the  development  of  an 
optimal  control  problem  to  the  application  of  optimal  control  to  this  research.  Chapter  II 
provides  an  introduction  of  optimal  control  and  the  process  that  is  used  to  solve  an 
optimal  control  problem.  Chapter  III  introduces  the  launch  problem  that  is  to  be  solved 
by  this  thesis  and  also  provides  the  hand  calculations  that  set  up  the  boundary  value 
problem.  These  are  used  later  for  verification  and  validation  of  the  pseudospectral 
optimal  control  solution.  Chapter  IV  first  starts  off  with  a  validation  of  the  results  to 
demonstrate  that  an  optimal  solution  has  been  found.  The  chapter  then  displays  the  results 
for  visualization  of  the  trajectory  and  to  prove  that  the  results  obtained  were  optimal 
using  the  derived  equations  from  Chapter  III  and  propagation  of  the  controls.  Chapter  IV 
then  introduces  the  uncertainty  analysis  that  was  performed  to  test  for  variations  in 
environmental  parameters.  Chapter  VI  gives  some  conclusions  and  suggests  some  ideas 
for  future  work. 
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II.  OPTIMAL  CONTROL  THEORY 


A,  INTRODUCTION 

The  popularity  for  using  optimal  control  is  based  on  three  main  reasons.  The  first 
reason  is  that  there  is  a  cost  function  associated  with  the  problem  that  can  be  minimized. 
The  types  of  minimized  cost  can  be  time,  fuel,  effort,  or  any  other  performance  objective. 
As  stated  before,  the  objective  of  this  thesis  research  is  to  maximize  final  velocity.  The 
next  benefit  to  optimal  control  is  the  use  of  dynamics  equations.  The  dynamics  equations 
allow  the  user  to  more  accurately  model  a  trajectory  that  the  system  can  fly.  As  compared 
to  kinematics  only,  this  allows  for  the  prediction  of  what  the  system  will  do  to  a  high 
degree  of  accuracy.  Lastly,  optimal  control  provides  the  ability  to  apply  constraints  to  the 
system  to  be  able  to  control  the  behavior.  The  constraints  can  be  in  the  form  of  time, 
states,  controls,  and  boundary  conditions. 

The  process  that  is  used  to  solve  optimal  control  problems  involves  first 
constructing  the  Hamiltonian.  The  Hamiltonian  is  key  to  deriving  the  Hamiltonian 
minimization,  the  costate  dynamics  (adjoint)  equations,  and  the  transversality  condition. 
In  the  following  analysis,  x  is  the  state  variable  and  u  is  the  control  variable.  The  costate 

vector  used  in  this  analysis  is  the  Lagrange  multiplier  function  and  is  annotated  as  yi(t) 
[10].  The  addition  of  this  function  will  be  further  described  in  the  next  section. 
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B.  A  GENERIC  OPTIMAL  CONTROL  PROBLEM 

A  generic  optimal  control  problem  (OCP)  is  given  as: 

x(t)  =  f[x[t),u[t)^ 

=  (1) 
t  =f 

e{x(l,)]  =  0 

In  (1),  /[x(.),m(.)]  is  the  cost  function  that  should  be  minimized.  The  cost  function  is 

composed  of  the  endpoint  cost,£’ ^  and  the  miming  cost,  j  ^  F(x(t),M(t)yt . 

The  endpoint  cost  is  associated  with  the  final  time  of  the  simulation.  Example  endpoint 
costs  can  be  final  time,  remaining  fuel,  terminal  velocity,  etc.  The  mnning  cost  is  cost 
accumulated  during  the  entire  flight  time.  An  example  of  mnning  cost  is  control  effort. 

The  dynamics  portion  is  defined  byi  =  the  initial  condition  is  defined  as 

x“ ,  and  the  start  and  end  times  are  defined  by  t  °  and  t  ^  .  Lastly,  any  endpoint  constraints 
are  contained  in  the  equation  ))  =  0  •  An  example  of  an  endpoint  constraint  can  be 

the  conditions  to  maintain  a  specific  orbit. 


Figure  1.  Optimal  maneuver,  from  [11] 


The  application  of  the  above  equations  is  seen  in  Figure  1  as  the  object  starts  at 


some  initial  condition  and  maneuvers  along  a  trajectory  to  a  desired  end  condition  along  a 
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given  dynamic  constraint.  The  problem  starts  at  and  the  spacecraft  maneuvers  itself 
to  endpoint,  which  satisfies  the  constraint  of  )  =  0 


C.  SOLVING  AN  OPTIMAL  CONTROL  PROBLEM 

In  the  1950s,  finding  solutions  to  the  standard  problem  formulation  given  in  (1) 
was  causing  problems  for  Soviet  engineers  who  knew  that  the  problems  being 
encountered  were  from  the  math  and  not  the  engineering.  This  led  the  Russian  military  to 
approach  an  individual  by  the  name  of  Lev  Pontryagin  to  help  solve  this  problem.  While 
creating  a  general  problem  of  optimal  control,  Pontryagin  realized  that  constraints  on  the 
control  need  to  be  included  and  special  attention  needs  to  be  given  to  optimization.  This 
gave  birth  to  the  present  form  of  optimal  control  theory  [12]. 


It  is  the  minimization  of  the  Hamiltonian  that  needs  to  be  given  the  proper 
attention.  When  the  Hamiltonian  is  minimized,  the  endpoints  of  the  control  constraint 
need  to  be  evaluated  to  determine  true  minimum.  This  process  is  described  as: 

rmnH[A,x,u) 


<u<u^ 


(2) 


Equation  (2)  is  the  Hamiltonian  minimization  condition  (HMC).  Pontryagin 
proved  that  the  minimized  Hamiltonian  is  always  constant  as  a  function  of  time  with  the 
value  zero  for  problems  independent  of  time,  -1  for  a  minimum  time  problem  [12]. 


Solving  Pontryagin’s  problem  can  be  decomposed  into  four  steps.  When  solving 
the  problem,  as  stated  above,  the  first  step  is  to  solve  for  the  Hamiltonian.  The 
Hamiltonian  is  a  function  of  the  running  cost  and  the  product  of  the  costate  vector  with 
the  dynamics  equations: 

H[A.,x,u)\=  F[x,u)  + f[x,u)  (3) 

The  next  step  is  to  perform  the  Hamiltonian  minimization.  This  involves  taking 
the  partial  derivative  of  the  Hamiltonian  with  respect  to  each  control  variable  and  setting 
the  derivative  equal  to  zero. 


The  goal  of  this  step  is  to  be  able  to  remove  the  dependence  of  the  control 
variable,  u,  from  the  Hamiltonian  equation.  If  u  does  not  appear  explicitly,  the  partial 
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derivative  is  interpreted  as  a  switching  function.  The  switching  function  is  describes  how 
u  switches  from  ui  to  Uu  (the  control  bounds)  throughout  the  maneuver.  The  next  step  is  to 
derive  the  adjoint  equations.  This  step  forms  the  dynamics  of  the  costate  variables,  as  a 
function  of  time,  for  each  state  variable.  The  adjoint  equation  is  needed  because  the 
minimized  Hamiltonian  is  a  function  of  X{t)  Therefore,  the  costate  needs  to  be  solved. 

This  is  done  by  taking  the  negative  of  the  partial  derivative  of  the  Hamiltonian  with 
respect  to  each  of  the  state  variable.  This  gives  the  adjoint  equation. 

The  last  step  is  to  apply  the  transversality  condition.  This  involves  taking  the 
partial  of  the  endpoint  Lagrangian  with  respect  to  the  state  at  the  terminal  time. 

(5) 

In  (5)  E  is  the  Endpoint  Lagrangian. 

E  :=  E  [x{tf  ))+  v'^e[x{tf  (6) 

From  here,  the  classical  approach  is  to  form  a  boundary  value  problem  (BVP) 
using  the  dynamics  and  adjoint  equations  together  with  the  boundary  conditions  and  the 
transversality  condition.  A  common  approach  to  solve  the  BVP  is  to  use  a  shooting 
method. 

D,  EXAMPLE  PROBLEM 

To  illustrate  the  application  of  the  idea  in  the  previous  section,  a  simple  example 
problem  will  be  studied.  The  example  that  will  be  solved  in  the  section  is  a  1-D  linear 
quadratic  problem  [11]. 
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J 


=^1, /'(')* 


x  =  x  +  u 

x„  =  l 

?o  =  0 

V  =  1 

X,  =  0 


(V) 


Examining  the  cost  function,  it  can  be  seen  that  e  {^x  {t^ ))  =  0  and  F  (x,  m  )  =  — • 
There  is  only  one  dynamics  equation  and  that  is  x  =  f[x,u^  =  x  +  u.  From  here,  the 


Hamiltonian  can  now  be  derived. 


//(T,x,m)  =  F(x,  m)  +  T^/(x,  m)  =  +  T(x  +  m)  (8) 

Now,  the  Hamiltonian  minimization  is  accomplished  by  taking  the  partial 
derivative  of //with  respect  to  u. 

^  =  M  +  A  =  0  (9) 

du 

This  allows  for  u  to  be  solved  for  in  terms  of  X,  which  gives  u  =  -A .  Since  w  is  a 
function  of  A(t)  ,  the  adjoint  equation  will  be  needed  to  solve  for  the  costate  history.  The 
adjoint  equation  is  the  next  step  in  solving  the  problem.  This  involves  taking  the  partial 
derivative  of  the  Hamiltonian  with  respect  to  x. 


dx 

A  =  -A 


(10) 


From  (10),  it  can  be  seen  that  the  solution  to  the  adjoint  equation  X  is  an 
exponential.  Now,  the  last  part  is  to  apply  the  transversality  condition.  The  first  part  of 
this  involves  solving  for  the  endpoint  Fagrangian. 

E  :=  E  [x{t f  e{^x{t ^  =  ux f  (11) 

Once  the  endpoint  Fagrangian  is  constructed,  to  obtain  the  transversality 
condition,  the  partial  derivative  of  the  endpoint  Fagrangian  is  taken  with  respect  to  the 
endpoint  condition. 


9 


=  V 


(12) 


dE 

dXf 


The  result  in  (12)  shows  that  the  value  of  the  costate  is  an  unknown,  u,  at  t^.. 

Thus  no  new  information  results  from  this  step.  The  transversality  condition  is  not 
necessary,  in  this  case,  because  only  two  boundary  conditions  are  needed  and  these  are 
provided  by  the  given  problem.  After  these  four  steps  are  completed,  the  following 
boundary  value  problem  can  be  constructed. 

X  —  X  —  A. 

A- -A 


Xq  =1 

x^  =  0 


(13) 


From  here,  the  problem  can  be  solved  numerically  to  obtain  A(t)  and  hence  u(t), 
which  is  desired.  This  process  is  not  necessarily  easy  to  accomplish.  Due  to  the  instability 
of  the  Hamiltonian  system,  the  integrated  equation  can  “blow  up”  even  in  the  face  of  a 
very  accurate  guess  for  the  unknown  initial  values  [13].  This  is  where  the  MATLAB  tool 
DIDO  can  make  life  a  lot  easier.  Once  the  cost  function,  dynamics  equations,  constraints, 
and  events  are  programmed  into  DIDO,  the  algorithm  will  solve  for  the  states,  controls, 
Hamiltonian,  and  costates  as  a  function  of  time,  without  the  need  to  construct  the  BVP. 
This  is  much  easier  than  building  for  the  BVP  and  using  a  shooting  algorithm  to  solve  the 
problem.  It  takes  away  the  need  to  build  an  algorithm  that  converges  on  a  solution 
without  an  accurate  initial  guess.  With  the  optimal  control  trajectories,  they  can  be 
propagated  to  solve  for  the  states  via  an  ordinary  differential  equation  (ODE)  for 
verification  and  validation  purposes. 


E, 


SUMMARY 


This  chapter  explained  why  optimal  control  is  widely  used  based  on  its  many 
advantages.  It  then  went  on  to  set  up  a  generic  OOP  that  was  to  be  solved  using  the 
method  defined  in  this  chapter.  After  the  process  for  solving  an  OCP  was  defined,  an 
example  problem  was  introduced  to  further  illustrate  the  procedure.  The  next  chapter  will 
define  the  launch  vehicle  problem  that  is  to  be  addressed  by  this  thesis  research. 
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III.  LAUNCH  PROBLEM  FORMULATION 


A,  INTRODUCTION 

In  this  chapter,  the  launeh  vehiele  problem  is  defined  along  with  the  desired  goal. 
Here  Pontryagin’s  prineiple  is  used  to  set  up  the  BVP  but  it  is  not  eompletely  solved  in 
this  ehapter.  The  BVP  provides  information  that  ean  be  ehecked  to  verify  that  an  optimal 
solution  has  been  found.  In  the  next  ehapter,  the  problem  is  solved  using  DIDO. 

B.  THE  LAUNCH  PROBLEM 

This  seetion  will  identify  the  variables  and  parameters  that  will  be  used  to 
eonstruet  the  launeh  OCP  starting  with  states,  eontrols,  eost,  and  lastly  the  dynamics 
equations.  The  first  part  to  building  this  problem  is  to  define  the  state  veetor  and  the 
eontrol  veetor.  The  state  veetor  for  the  problem  ineludes  the  Cartesian  positions, 
veloeities,  and  thrust  direetion  eosines.  They  are  shown  in  Equation  14. 


The  positions  (x,  y,  and  z)  are  in  units  of  kilometers  (km),  the  veloeities  (vx,  Vy,,  and  v^) 
are  in  kilometers/seeond  (km/s),  and  the  thrust  direetion  eosines  (cx,  Cy,  and  Cz)  are  unit 
less  as  this  unit  veetor  that  simply  provides  the  direetion  of  the  eonstant  thrust.  The 
eontrols  are  the  rates  of  ehange  of  the  unit  thrust  veetor  and  the  Euelidean  distanee  of  the 
launch  vehicle  from  the  origin  of  a  referenee  frame.  The  radius  veetor  was  added  as  a 
eontrol  to  introduee  a  eonstraint  to  prevent  the  launeh  vehiele  from  entering  the  surface  of 
the  Earth,  i.e.  r>r^.  The  eonstraint  r  >  also  avoids  a  potential  situation  where  zero 
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would  be  in  the  denominator  of  the  dynamies  equation  (see  Equation  17).  The  eontrol 
vector  is: 


U  = 


^y 

r 


(15) 


The  goal  for  this  problem  is  to  maximize  the  final  velocity  of  the  launch  vehicle  and  this 
becomes  the  cost  function.  Final  velocity  was  chosen  to  be  the  endpoint  cost  due  to  the 
desire  to  achieve  maximum  velocity  in  the  quickest  manner  possible.  The  bounds  on  time 
that  were  used  were  a  starting  time  of  zero  and  a  final  time  based  on  how  long  it  took  to 
consume  all  first  stage  propellant.  This  allows  for  the  launch  vehicle  to  be  closer  to  final 
orbital  velocity  at  first  stage  burnout.  At  a  higher  first  stage  burnout  velocity,  less  thrust 
input  is  required  from  the  second  stage  to  achieve  the  desired  orbit.  When  performing  the 
analysis,  the  cost  function  is  the  variable  that  is  to  be  minimized.  In  order  to  maximize 
the  final  velocity,  the  cost  function  has  to  be  the  negative  of  the  final  velocity,  which  is 
the  same  as  minimizing  the  negative  of  final  velocity  seen  in  Equation  16.  The  square  of 
velocity  was  used  to  remove  the  need  to  have  a  square  root  in  the  equation.  This 
eliminates  the  potential  of  having  a  square  root  of  zero,  which  has  an  infinite  gradient. 

/[x(»), «(•)]  = -v/  (16) 

We  now  define  the  dynamics  equations  which  will  govern  this  problem.  These  equations 
are  formed  from  the  time  rate  of  change  of  the  state  variables.  The  dynamics  are  shown  in 
Equation  17. 
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In  (17),  T  is  the  thrust  of  the  launeh  vehicle  in  N  which  is  held  constant  during  the  flight 
time,  m  is  the  mass  of  the  vehicle  in  kg,  Isp  specific  impulse  of  the  vehicle  in  sec,  //  is  the 

3  2  3  2 

gravitational  constant  of  the  earth  in  kg  /s  ,  /?  is  the  atmospheric  density  in  kg/m  ,  v 

is  the  relative  velocity  of  the  vehicle  with  the  atmosphere  in  km/s,  S  is  the  surface  area  of 
the  vehicle  in  m  ,  and  Cd  is  the  coefficient  of  drag. 

A  path  constraint  was  added  to  the  problem  in  order  to  maintain  the  magnitude  of 
the  thrust  direction  cosines  equal  to  one  and  the  radius  from  the  states  x,  y,  and  z  is  equal 
to  the  control  radius.  These  constraints  were  necessary  to  ensure  that  the  defined 
magnitude  of  thrust  would  not  be  exceeded  and  launch  vehicle  flight  path  remained 
outside  the  radius  of  the  Earth.  Than  path  constraints  are  given  as: 


c/  +  cj+c,^-l  =  0 

2,2,2  2  A 

X  +  y  +z  -r  =0 


(18) 
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The  full  optimal  control  problem  is  now  given  as: 

j[x(.),w(.)]  =  -v/ 
x  =  v^ 

Z  =  V. 


T  2 

V.  =-c, — - 

m  r 


m 


T  n  2 

'>=^^>'-7^ - - 


1  2 


T  M  1 

=  — C. - - 

m  r 


P^\elJC, 


m 


m  =  ■ 


Ispg 


c  =w 

X  X 


Cv  =  Wv 

y  y 


Xg  =  cos{Lat)cos{Lon) 
Jo  =  cos(Tat)  sm(Ton) 
Zg  =  sm(Tat) 
c/+c7c,^-l  =  0 
x^  +  +z^  -r^  =0 

to  =0 


^/  = 


Mq  -  rUf 

m 


(19) 


In  (19),  Lat  and  Lon  are  the  latitude  and  longitude  of  the  launch  point  and  is  the  radius 
of  the  Earth. 


C.  DEVELOPING  THE  BOUNDARY  VALUE  PROBLEM 

As  described  in  the  previous  chapter,  the  first  step  in  setting  up  the  OCP  BVP 
involves  solving  for  the  Hamiltonian.  Because  the  cost  is  only  a  function  of  final  velocity, 

the  running  cost  is  zero  and  therefore  F  (^x,u^  is  zero.  The  Hamiltonian  is  now  just  a 
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function  of  the  individual  costates  and  time  rate  of  change  of  each  state  variable. 

A  'r  1  \ 

’'y'y  ' ''z’ z  '  '  "vx  I  ''x  3  r\  ’  rel ,x'~’ d 


V 


m 


r 


/I, 


vy 


S  3  ^  P^rel,y^^d 


f  • 


+  A.. 


jj.z  1 


ym  r  2  j 


(20) 


\m  '  r“  2 
L K  +  \ %  +  \ (L  +  <^^  +  <^z  “ l)  +  L‘>a,r 

Now  that  the  Hamiltonian  is  formed,  the  next  part  is  to  perform  the  Hamiltonian 
minimization.  The  partial  derivative  of  H  is  performed  with  respect  to  each  of  the 
controls.  The  resulting  relationships  are  given  in  (21). 
dH 


dw^ 

dH 

dw^, 

dH 

dw. 


.=  2  =0;=5, 


=  A  =0;=5, 


=  ;  =0;=5, 


(21) 


dH  3/1^ /iz 


dr 


-^K^,h,y  =  ^ 


As  seen  in  (21),  the  Hamiltonian  is  linear  in  the  thrust  direction.  Therefore,  in  accordance 
with  Pontryagin’s  principle  defined  in  the  previous  chapter,  the  partials  need  to  be 
interpreted  as  switching  functions  shown  by  5'j,5'2,  and5'3 .  The  adjoint  equations  are 
constructed  next  by  taking  the  partial  derivative  of  H  with  respect  to  the  state  vectors. 
The  relationships  are  given  by  (22). 
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5  _dH  .. 

^  3  ^'^path,r^ 

OX  r 

'^y  =  ^  =  ^-^Kath,ry 
oy  r 

,  _dH  .. 

~  ^  ~  3  ■^‘^pathy^ 

oz  r 

K  =^  =  -4 
A.  =^  =  -4 

a/f  -KT 

^liy  ^  ^path,u^x 

oc^  m 

dH  T 

1.  =  — = 


5c„ 


m 


•  _dH  _-KT 


(22) 


dc. 


m 


For  the  transversality  condition,  the  endpoint  Lagrangian  is  based  completely  on 
the  endpoint  cost  of  maximizing  final  velocity. 


(23) 


Now  the  partial  derivative  of  the  endpoint  Lagrangian  is  performed  for  each  of  the 
velocities  and  those  are  shown  in  Equation  24. 

dE 


^.,(4)  =  ^  =  -2''4'/) 
L  (4) =^=-24(4) 

4,('/)  =  ^  =  -24((,) 


(24) 


In  (24),  it  can  be  seen  that  the  velocity  costate  endpoint  is  related  to  the  final  velocity. 
This  can  be  useful  as  a  verification  and  validation  result.  Similarly,  in  (22)  the  velocity 
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adjoint  is  a  function  of  the  position  costate.  This  suggests  that  those  costates  may  vary 
linearly  whieh  is  also  useful  as  a  eheek  during  the  verifieation  and  validation. 

D.  SUMMARY 

This  chapter  defined  the  launch  problem  and  showed  how  Pontryagin’s  principle 
can  be  used  to  construct  the  BVP  for  the  launch  problem.  Once  the  BVP  is  constructed,  it 
would  be  a  very  ehallenging  proeess  to  obtain  a  solution  using  a  shooting  method  (e.g., 
POST).  The  results  obtained  here  will  be  used  to  verify  a  candidate  solution.  In  the  next 
ehapter,  DIDO  is  used  to  solve  this  problem  with  the  goal  of  reducing  computational 
time. 
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IV.  LAUNCH  TRAJECTORY  OPTIMIZATION 


A.  INTRODUCTION 

In  this  chapter,  an  optimal  control  solution  is  obtained  with  DIDO  and  verification 
and  validation  of  the  results  is  performed  to  indicate  optimality  of  the  solution.  Then,  the 
individual  results  will  be  displayed  to  analyze  and  illustrate  the  trends. 

B.  OPTIMIZATION  RESULTS 

The  parameters  that  were  used  in  this  problem  formulation  can  be  seen  in  Table  1. 
After  running  DIDO,  an  optimal  solution  for  16  nodes  was  found  indicating  the  problem 
was  correctly  posed.  To  confirm  the  results  from  the  output  of  DIDO  and  series  of  plots 
were  created  for  verification. 


Parameter 

Value/Range 

mo 

219676  kg 

mf 

6145  kg 

Isp 

397.45  sec 

T 

960000  N 

S 

2.17  m^ 

Cd 

0.15 

X,  y,  z 

-6800  to  6800  km 

Vx,  Vy,  Vz 

-5  to  5  km/s,  -5  to  5km/s,  0  to  5  km/s 

Cxj  Cy,  Cz 

-1  to  1,  -1  to  1,  0  to  1 

Wx,  Wy,  Wz 

-0.1  to  0.1 

r 

6378.1363  to  6800  km 

Lat 

28°N  31  min  26.61  sec 

Lon 

80°W  39  min  3.06  sec 

Table  1.  Model  parameters,  from  [14] 
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After  examining  Figures  2  and  3,  the  output  is  what  is  desired  in  that  there  is  a  parabolie 
inerease  in  the  launeh  vehiele’s  veloeity  whieh  obtains  a  final  value  of  4.73  km/s.  The 
magnitude  of  the  unit  thrust  veetor  is  also  eonstant  at  unity  as  required. 


Figure  2.  Position  trajeetories  for  maximum  final  veloeity 


time  (sec) 

Figure  3.  Veloeity  and  unit-thrust  veetors  for  maximum  final  velocity 
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After  reviewing  the  plots  of  the  eostates  from  Figure  4,  they  behave  as  expected 
based  on  the  adjoint  equations  that  were  derived  from  (22).  Specifically  examining  the 
adjoint  equations  for  the  velocity  costates,  the  derivative  of  the  individual  costate  is  the 
negative  of  the  position  costate. 


1  1 

'  1 

1  1  — _ 

1  1 

1  1 

1  1 

1  - — 

^vx 

— 1 

vy 

^vz 

1  1 

1  1 

■1.5 1 - 1 - 1 - 1 - 1 - 1 - 1 - - 
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Figure  4.  Costate  trajectories  for  maximum  final  velocity 


In  Figure  5,  looking  at  ,  it  starts  with  a  positive  slope  and  has  a  maximum 
around  45  seconds.  The  plot  of  -  starts  off  positive  and  crosses  the  y  axis  at  the  same 
time  X^  attains  a  maximum.  The  same  can  be  done  for  the  adjoint  equations  of  Vx  and  Vy. 
As  per  (24),  the  final  value  of  X^,^,  X^,  and  2,,^  are  required  to  be  )  =  -8.92km  /  5 , 

-2Vy{ty)  =  -1.08km  /  5  ,  and-2v^(t^)  =  -2.91km  /  s ,  respectively,  which  is  not  the  case 
after  examining  Figure  5.  This  inconsistency  comes  from  the  scaling  that  was  used  while 
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running  DIDO.  After  removing  the  veloeity  sealing  of  7.9054  km/s,  the  final  veloeity 
costate  values  agreed  with  the  transversality  conditions.  This  analysis  further  validates 
the  results  obtained  by  DIDO. 


Figure  5.  Verification  and  validation  of  costates 


Examining  Figure  6,  the  altitude  of  the  launch  vehicle  increases  as  the  velocity 
increases  and  at  burnout,  attains  an  altitude  of  33.59  km  above  the  surface  of  the  earth. 
The  rapid  rise  in  altitude  is  expected  as  the  launch  vehicle  is  at  a  constant  thrust  which  is 
creating  a  linear  rise  in  velocity  coupled  with  an  exponential  decrease  in  atmospheric 
density  reducing  aerodynamic  drag  encountered  by  the  launch  vehicle. 
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time  (sec) 


Figure  6.  Altitude  as  a  function  of  time 

Figure  7  represents  the  control  vectors  obtained  from  the  solution.  The  radius  was 
removed  from  the  plot  to  be  able  to  view  the  thrust  rate  of  change  more  accurately.  The 
thrust  rate  peaked  at  0.1  at  the  very  beginning  which  is  the  maximum  defined  by  Table  1 
but  tapered  off  as  time  increased  and  eventually  converged  to  zero.  Since  this  is  not  a 
minimum  time  problem,  the  Hamiltonian  constant  is  required  to  be  zero.  Looking  at 
Figure  8  the  Hamiltonian  is  approximately  zero.  This  is  indicative  of  an  optimal  solution. 
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Figure  7.  Control  vectors 


Figure  8.  Hamiltonian  evolution  as  a  function  of  time 
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c. 


VERIFICATION  AND  VALIDATION 


In  order  to  verify  the  results  that  were  obtained  from  DIDO,  a  simulation  was  run 
using  the  ode45  solver  in  MATLAB.  The  eontrol  veetor  whieh  was  obtained  via  DIDO  is 
used  to  propagate  the  solution  to  verify  the  results.  This  will  eonfirm  whether  or  not  the 
solution  obtained  via  DIDO  is  feasible  for  implementation.  When  the  two  results  are 
plotted  against  eaeh  other,  it  can  be  seen  if  the  trends  are  the  same  or  if  there  are  large 
disparities  in  the  data.  If  the  two  results  are  the  same,  it  shows  that  the  optimal  solution 
obtained  from  DIDO  is  a  valid  one.  If  they  are  different,  then  the  DIDO  solution  may  not 
be  accurate  enough  and  the  problem  set  up  may  need  to  be  reevaluated  and  solved  with 
larger  number  of  nodes.  Figure  9  shows  the  plots  of  the  DIDO  solution  compared  with 
the  propagated  solution 
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Figure  9.  Verification  and  validation  of  DIDO  solution 
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In  Figure  9,  it  can  be  seen  that  the  two  solutions  are  nearly  identical.  This  shows  that  the 
solution  obtained  from  DIDO  is  indeed  a  feasible  one. 

D.  FURTHER  ANALYSIS 

After  the  optimal  results  were  obtained  from  DIDO,  there  was  some  further 
analysis  done  to  assist  with  better  visualization  of  the  trajectory.  The  first  analysis  was  a 
coordinate  transformation  as  seen  in  Figure  10.  The  simulation  is  best  solved  in  ECI 
coordinates  since  most  launches  target  a  specific  orbit  and  it  is  best  to  maintain  the  state 
vectors  in  a  coordinate  system  that  is  centered  intertially  in  the  Earth.  To  better  visualize 
the  trajectory  of  the  launch  vehicle  a  rotation  matrix  was  applied  to  the  coordinates  to 
transform  them  from  Earth  centered  inertial  to  a  north-west-up  frame  shown  in  Eigure  10. 
This  was  done  by  rotating  about  the  z-axis  to  position  the  y-axis  at  the  longitude  of  the 
launch  point.  Next,  the  coordinate  system  was  rotated  about  the  new  x-axis  to  point  the  z- 
axis  to  the  latitude  of  the  launch  point.  Eastly,  the  coordinate  system  was  rotated  about 
the  new  z-axis  to  point  the  x-axis  in  the  north  direction.  The  resulting  coordinate  system 
now  has  the  x-axis  pointing  north,  y-axis  point  west,  and  the  z-axis  as  the  zenith 
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Figure  10.  Coordinate  transformation  from  ECI  to  NWU 


The  transformation  matrices  used  are  given  in  (25)  where  Lon  is  the  longitude  of 
the  launch  site  and  Lat  is  the  latitude. 
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After  substituting  the  Lon  and  Lat  for  KSC,  we  obtain  a  final  rotation  matrix. 

■-0.0776  0.4712  0.8786' 

C  =  C3C2  Cl  =  -0.9867  -0.1625  0  (26) 

0.1427  -0.8669  0.4775 

Figures  11  and  12  show  the  launeh  vehiele’s  trajeetory  after  performing  the  eoordinate 
transformation  with  the  origin  of  the  new  eoordinate  system  being  the  launeh  point.  The 
results  are  exactly  as  expected  in  that  initially  the  thrust  vector  is  vertical  direction  and 
beginning  to  turn  over  into  the  direction  of  flight.  That  is  consistent  with  the  constraint 
that  was  applied  to  ensure  that  the  initial  thrust  vector  was  aligned  with  launch  point’s 
radius.  The  velocity  is  at  all  times  tangential  to  the  launch  vehicle’s  trajectory. 
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Figure  1 1 .  Launch  vehicle  trajectory  with  thrust  vectors 
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Figure  12.  Launch  vehicle  trajectory  with  velocity  vectors 
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The  next  analysis  performed  was  to  plot  the  rotated  trajectory  on  a  Google  Earth 
plot  to  show  how  the  trajectory  performed  when  plotted  in  reference  to  the  land  mass. 
Figures  13  and  14  provide  a  perspective  of  the  launch  vehicle  as  it  leaves  the  launch  point 
and  travels  in  a  north  easterly  direction.  This  is  consistent  with  launches  that  are  currently 
done  at  the  KSC  in  that  once  the  launch  vehicle  leaves  the  launch  pad,  it  heads  in  a  north 
easterly  direction  to  ensure  a  safe  area  to  drop  the  boosters.  Figure  15  shows  a  plot  of  the 
STS-135  launch  [15].  It  can  be  clearly  seen  that  the  trajectory  developed  using  DIDO  is 
very  similar  to  the  trajectories  used  by  NASA  for  the  space  shuttle  (shown  by  the 
trajectory  given  in  Figure  15). 
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Figure  14.  Google  earth  2D  view  of  the  launch  trajectory 


Figure  15.  Google  earth  2D  view  of  STS-135  trajectory,  from  [15] 


E.  PROBLEMS  ENCOUNTERED 

Throughout  the  process  of  creating  a  solution  to  this  problem,  there  were  many 
hurdles  that  had  to  be  overcome  in  order  to  produce  viable  results.  The  first  challenge 
was  establishing  a  type  of  coordinate  system  that  was  to  be  used  in  order  to  model  the 
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launch  problem.  First,  the  problem  formulation  was  modeled  in  Cartesian  eoordinates. 
This  ereated  problems  in  effeetively  being  able  to  maintain  the  trajeetory  from  eolliding 
with  the  surfaee  of  the  Earth.  The  problem  formulation  was  then  shifted  to  polar 
eoordinates  whieh  allowed  for  maintaining  the  radius  of  the  trajeetory  outside  the  surfaee 
of  the  earth.  The  solution  was  then  obtained  in  drag  free  environment  but  it  was  still 
desired  to  keep  the  problem  in  Cartesian  eoordinates.  The  problem  was  then  shifted  baek 
to  Cartesian  and  an  optimal  solution  was  finally  obtained  for  a  drag  free  environment. 
Onee  drag  was  introdueed  it  was  beeoming  impossible  to  keep  the  launeh  vehiele  from 
eolliding  with  the  surfaee  of  the  earth.  After  extensive  isolation  of  the  components  to  the 
dynamies  equation,  it  was  determined  that  the  equation  in  whieh  density  was  being 
ealculated  was  the  souree  of  the  problem.  Once  that  was  isolated,  an  optimal  solution  was 
found  which  lead  to  ereating  the  method  of  solving  for  density  that  will  be  mentioned  in 
the  next  ehapter.  Overeoming  these  ehallenges  emphasizes  that  proper  problem 
formulation  is  eritieal  to  suoeessfully  solving  the  problem. 

F.  SUMMARY 

In  this  ehapter,  the  optimal  eontrol  solution  for  launeh  was  presented  and 
evaluated  for  feasibility.  The  validation  eheeks  were  aeeomplished  using  the  derived 
equations  from  the  previous  ehapter.  Onee  those  eheeks  were  eomplete,  the  eontrols 
obtained  via  DIDO  were  propagated  to  obtain  a  new  set  of  state  variables.  The  two  sets  of 
states  were  plotted  against  eaeh  other  to  establish  feasibility  of  the  solution.  The  two 
results  were  nearly  identieal  proving  the  solution  was  feasible.  Next,  the  trajeetory  was 
transformed  to  another  coordinate  frame  for  visual  reference  and  plotted  using  a  Google 
Earth  to  evaluate  how  the  trajeetory  performs  with  land  mass  visible.  Onee  again  that  was 
ehecked  against  trajeetories  flow  by  NASA  for  the  spaee  shuttle.  The  present  solution  is 
very  similar  to  existing  trajeetories. 
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V.  MONTE  CARLO  SIMULATION 


A,  INTRODUCTION 

One  way  to  assess  the  impact  of  uncertainties  in  a  dynamic  system  is  by  the  use  of 
a  Monte  Carlo  analysis.  This  is  done  by  drawing  from  a  large  pool  of  random  samples 
and  observing  their  behavior  [16].  This  chapter  starts  off  by  describing  the  density  model 
used  to  model  atmospheric  density  as  a  function  of  temperature  offset  and  altitude.  Next 
the  Monte  Carlo  simulations  that  were  performed  to  assess  the  effects  of  uncertainties  in 
launch  environmental  are  described.  Three  Monte  Carlo  simulation  studies  were 
performed.  The  first  was  a  temperature  only  simulation,  the  second  was  a  wind  only 
simulation,  and  the  third  was  a  combination  of  both  temperature  and  wind.  Temperature 
and  wind  were  chosen  based  on  the  possibility  of  these  effects  having  the  largest 
influence  on  the  launch  vehicle. 


B,  DENSITY  MODELING 


One  model  predicting  atmospheric  density  as  a  function  of  altitude,  r,  uses  an 
exponential  form  [14]  . 


Pir)  =  p^QW 


Tq  —  r  ^ 
H 


(27) 


In  (27),  po  is  the  atmospheric  density  at  sea  level  and  H  is  the  scale  height  parameter.  For 
the  density  analysis  in  this  thesis,  various  data  points  were  taken  from  the  1976  Standard 
Atmosphere  in  order  to  calculate  density  as  a  function  of  altitude  {alt)  and  temperature 
[17].  Data  points  were  taken  for  temperature  offsets  {TO)  from  -30°C  to  30°C  in  10°C 
increments.  Once  those  points  were  obtained  they  were  plotted  on  an  Excel  scatter  plot 
and  a  sixth  order  polynomial  was  used  to  create  a  curve  fit  of  density  versus  altitude  for 
each  TO  as  seen  in  Figure  16.  The  format  of  the  equation  for  modeling  density  is  given  as 
p  -  aplf  +  GjaW  +  a^alt'^  +  a^alt^  +  a^alt^  +  a^alt  +  (28) 


Table  2  shows  the  coefficients  at  each  power  of  alt  for  the  given  TO. 
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Figure  16.  Air  density  as  a  function  of  altitude  and  temperature 


TO  (°C) 

ai 

a2 

a3 

a4 

as 

ae 

a? 

-30 

9.482e-‘' 

-9.176e-" 

3.329e-^ 

-5.958e-^ 

8.025e^ 

-0.1355 

1.370 

-20 

8.699e-‘' 

-8.380e'" 

3.028e'^ 

-5.444e'^ 

7.587e^ 

-0.1312 

1.316 

-10 

8.007e-" 

-7.680e-' 

2.765e'" 

-4.996e-^ 

7.199e' 

-0.1271 

1.271 

0 

7.398e-‘' 

-7.066e-' 

2.536e-^ 

-4.607e-^ 

6.855e^ 

-0.1233 

1.227 

10 

6.855e-‘' 

-6.523e-^ 

2.333e-^ 

-4.266e-'* 

6.546e^ 

-0.1197 

1.186 

20 

6.373e-" 

-6.042e-' 

2.155e-" 

-3.965e-^ 

6.269e' 

-0.1163 

1.147 

30 

5.944e-‘' 

-5.615e-" 

1.998e'^ 

-3.700e-^ 

6.019e^ 

-0.1131 

1.111 

Table  2.  Coefficients  from  density  curve  fits 


Each  of  the  columns  of  coefficients  given  in  Table  2  were  also  fit  to  a  linear  curve 
to  be  able  to  be  able  to  compute  each  coefficient  as  a  function  of  TO.  This  allows  a  single 
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equation  for  density  to  be  developed  as  a  function  of  alt  and  TO.  Figures  17  through  23 
show  the  curve  fits  for  the  given  values  of  TO,  for  each  of  the  coefficients  in  Table  2. 


X  10'® 


Figure  17.  Curve  fit  for  a.\ 


Figure  18.  Curve  fit  for  aa 
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X  10'® 


Figure  19.  Curve  fit  for  a3 
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Figure  20.  Curve  fit  for  a4 
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Figure  2 1 .  Curve  fit  for  as 


Figure  22.  Curve  fit  for  ae 
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Figure  23.  Curve  fit  for  a? 


Using  the  results  from  Figures  17  through  23,  a  single  equation  ean  be  developed 
for  ealeulating  density: 

p{alt,TO)  =  (-5.864e-"  *TO  +  7.5376''’ )  *  alt^  +(5.899e''  *TO-  7.2126'' )  *  alt^  + 
(-2.2046''  *ro  +  2.5926'')*a/t^  +(3.7366'®  *70-4.7056'") +  (29) 

(-3.3246'®  *  70  + 6.9296'®)  *a/t'  +(3.7296'"  *70  -  .1237)  *  a/t  -  4.3076'®  *  70  +  1.233 

Figure  24  shows  the  error  of  the  single  equation  (29)  and  exponential  density  (27)  with 
the  densities  obtained  from  the  1976  standard  atmosphere. 
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Figure  24.  Error  comparison  of  single  equation  and  exponential  for  TO=0°C 


The  error  corresponding  to  the  model  (29)  maintains  a  relatively  constant  value, 
close  to  zero  whereas  the  exponential  model  has  a  rather  high  error  during  the  first  two 
thirds  of  the  flight  regime.  Using  the  single  equation  that  was  produced  by  this  thesis,  the 
air  density  calculated  will  be  more  accurate  and  provide  a  better  estimate  of  atmospheric 
drag  during  flight.  Accordingly,  this  was  the  model  used  for  obtaining  the  optimal  launch 
trajectory  discussed  in  the  last  chapter. 

C.  TEMPERATURE  VARIATION 

After  the  solution  from  DIDO  was  validated  (see  Chapter  IV),  the  next  step  is  to 
add  a  certain  degree  of  variation  into  the  dynamics  so  the  effects  can  be  analyzed.  The 
first  simulation  analyzed  the  effects  of  temperature  variation.  A  Monte  Carlo  simulation 
was  run  for  1000  different  points  using  a  normal  temperature  distribution  to  introduce  the 
uncertainty.  The  1-a  variation  used  in  the  temperature  simulation  was  30°C.  Varying  the 
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temperature  changes  the  air  density  the  launch  vehicle  encounters,  which  changes  the 
amount  of  drag  felt  on  the  launch  vehicle.  The  equation  that  was  used  to  model 
temperature  variation  is  given  by  (30). 


TO  =  TO  +  (j,  *n 

norm  temp 


(30) 


In  (30),  TO  is  the  random  temperature  used  in  the  density  calculation,  TO  norm  = 
0°C  is  the  temperature  offset  based  on  current  conditions,  atemp  =  30°C  is  the  temperature 
variation,  and  n  is  a  random  number  produced  from  a  normal  distribution  with  ju  =  0  and 
fj  =  1 .  Figure  25  shows  the  resulting  trajectories  from  the  Monte  Carlo  simulation  while 
Figure  26  shows  the  endpoints  of  the  trajectories  in  the  north-up  plane.  The  variation  in 
the  north-west  plane  is  shown  in  Figure  27.  Referring  to  Figures  25  through  27,  it  appears 
that  the  launch  trajectory  is  quite  insensitive  to  large  variation  in  temperature  and  the 
desired  trajectory  can  still  be  achieved. 
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Figure  25.  Monte  Carlo  simulation  for  temperature  uncertainty 
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Figure  26.  Plot  of  endpoints  from  Monte  Carlo  for  temperature  in  North-Up 


Figure  27.  Plot  of  endpoints  from  Monte  Carlo  for  temperature  in  North-West 


A  plot  of  the  distribution  of  temperature  was  ereated  to  verify  that  a  eorreet 
distribution  was  being  used  during  the  Monte  Carlo  simulation.  That  distribution  is  seen 
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in  Figure  28.  The  distribution  looks  as  expected  with  a  mean  value  of  zero  and  a  standard 
deviation  of  30°C. 


Figure  28.  Temperature  offsets  for  Monte  Carlo  simulation 


D.  WIND  VARIATION 

The  next  part  of  the  Monte  Carlo  simulation  involved  introducing  a  wind 
variation.  The  expected  wind  patterns  were  obtained  from  the  NOAA  Earth  Systems 
Research  Laboratory  (ERSE)  for  the  periods  of  January  to  December  at  a  level  of  300  mb 
which  is  equivalent  to  30,000  ft.  [18].  Eigures  29  and  30  are  the  graphics  obtained  from 
ERSE  website  [18]. 
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Figure  29.  Wind  directions  over  North  America  at  the  surface,  from  [18] 


NCEP/NCAR  Rannalysia 


Figure  30.  Wind  directions  over  North  America  at  30,000  ft.,  from  [18] 

Figures  29  and  30  show  that  the  prevailing  winds  have  a  general  westerly 
direction  at  30,000  ft.  and  a  negligible  wind  component  at  the  surface.  For  the  Monte 
Carlo  simulation,  a  normal  distribution  was  used  with  a  1-a  variation  in  wind  magnitude 

of  wind  =17^^.  Once  the  winds  were  broken  down  into  components,  they  were 

multiplied  by  the  Gaussian  random  number,  n,  with  ju  =  0  and  cr  =  1  then  multiplied  by 
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the  inverse  of  the  transformation  matrix  given  by  (25).  This  was  used  to  convert  the  wind 
vector  to  ECI  to  be  used  in  the  dynamics  equations. 


windj^ 

wind^ 

- 

windjj 

0 

wind^ 

windf,, 

wind,. 

=  C" 

wind^ 

wind^ 

wind^j 

(31) 


The  X,  y,  and  z  components  of  the  wind  variation  were  the  added  to  the  relative  wind 
velocity  in  the  dynamics  equations  given  by  (32). 

Kela,ve  +  (32) 

In  (32),  Q  X  r  is  the  rotation  velocity  of  the  atmosphere  based  on  Earth’s  angular  rotation 
vector  (33)  and  the  current  value  of  r. 


Q  = 


0 

0 

7.2921 158553e^' 


{rad  /  s) 


(33) 


Eigures  31  through  32  are  the  resulting  plots  of  the  trajectories  from  the  Monte  Carlo 
simulation  and  the  trajectory  endpoints.  Eigure  31  is  the  full  trajectories  from  the  Monte 
Carlo  while  Eigures  32  and  33  are  the  trajectory  endpoints  in  north-up  and  north-west, 
respectively. 
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<-  West  (km)  North  ->  (km) 

Figure  3 1 .  Monte  Carlo  simulation  for  wind  uncertainty 


Figure  32.  Plot  of  endpoints  from  Monte  Carlo  for  wind  in  north-up 
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Figure  33.  Plot  of  endpoints  from  Monte  Carlo  for  wind  in  north-west 

The  variance  of  the  wind  in  the  north-up  plan  is  consistent  with  the  variance  in 
temperature  in  the  same  plane.  Wind  has  a  much  smaller  variance  in  the  West  direction 
which  could  suggest  that  the  wind  variation  has  a  smaller  effect  on  the  launch  trajectory 
than  temperature. 

E.  TEMPERATURE  AND  WIND  VARIATION 

The  last  part  of  the  Monte  Carlo  simulation  involved  testing  variations  in  both 
temperature  and  wind  described  by  (30)  and  (32).  Figures  34  through  36  are  the 
trajectories  that  were  obtained  and  the  trajectory  endpoints.  Figure  34  is  the  full  trajectory 
obtained  from  the  Monte  Carlo  while  Figures  35  and  36  are  the  trajectory  endpoints  in 
north-up  and  north-west,  respectively. 
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<-  West  (km)  North  ->  (km) 

Figure  34.  Monte  Carlo  simulation  for  temperature  and  wind  uneertainty 


Figure  35.  Plot  of  endpoints  for  Monte  Carlo  of  temperature  and  wind  in  north- 

up 
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Figure  36.  Plot  of  endpoints  for  Monte  Carlo  of  temperature  and  wind  in  north¬ 
west 

The  trajeetories  and  trajectory  endpoints  are  consistent  with  the  two  previous  Monte 
Carlo  simulations.  Tables  3  and  4  show  the  statistical  data  from  the  three  simulations. 
The  standard  deviation  of  the  position  endpoints  and  the  mean  final  velocity  are  the  two 
most  important  variables  that  need  to  be  evaluated. 


Simulation 

Mn  (km) 

Mw  (km) 

Mu  (km) 

<j^  (km) 

<Jw  (km) 

<j^  (km) 

Temp 

53.934 

-87.895 

33.621 

0.390 

0.402 

0.247 

Wind 

53.945 

-87.908 

33.630 

0.123 

0.011 

0.098 

Both 

53.949 

-87.898 

33.635 

0.269 

0.390 

0.150 

Table  3.  Endpoint  position  means  and  standard  deviations 
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Simulation 

(km/s) 

My„.  (km/s) 

My,  (km/s) 

f7„^(km/s) 

(km/s) 

(^y,  (km/s) 

Temp 

2.151 

-3.710 

1.116 

0.015 

0.020 

0.011 

Wind 

2.151 

-3.710 

1.117 

0.0035 

1.5870"’ 

0.002 

Both 

2.152 

-3.710 

1.117 

0.012 

0.020 

0.008 

Table  4.  Endpoint  velocity  means  and  standard  deviations 


From  Table  3,  it  can  be  seen  that  neither  of  the  uncertainties  considered  had  a 
significant  effect  on  the  mean  final  position,  but  the  temperature  uncertainty  gave  the 
largest  spread  of  endpoints.  This  suggests  that  more  emphasis  needs  to  be  placed  on 
predicting  temperature  than  predicting  wind  patterns.  In  Table  4,  neither  of  the 
uncertainties  had  any  appreciable  effect  on  the  mean  final  velocity  of  the  launch  vehicle. 
Similar  to  before,  the  wind  uncertainty  had  much  smaller  effect  than  temperature.  This 
further  confirms  the  conclusion  that  temperature  has  a  larger  influence  on  the  launch 
trajectory  than  wind. 

F.  SUMMARY 

In  this  chapter,  a  model  was  developed  to  more  accurately  predict  atmospheric 
density.  This  lead  to  the  formulation  of  a  single  equation  that  can  predict  atmospheric 
density  as  a  function  of  temperature  offset  and  altitude.  Next,  Monte  Carlo  simulations 
were  performed  to  assess  the  effects  of  uncertainties  in  the  model.  The  uncertain 
variables  chosen  were  temperature  and  wind  based  on  the  fact  that  these  quantities  have 
the  potential  to  produce  the  greatest  effects  on  the  launch  trajectory.  The  results  from  the 
simulations  were  plotted  as  trajectories  and  their  respective  endpoints.  The  statistical 
analysis  showed  that  temperature  had  the  largest  effect  on  the  launch  trajectory  of  the 
three  Monte  Carlo  simulations  that  were  run.  Based  on  the  results  from  this  chapter,  the 
nominal  solution  may  be  sufficient  for  control  of  the  first  stage. 
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VI.  CONCLUSIONS  AND  FUTURE  WORK 


A,  CONCLUSION 

The  stated  goals  for  this  thesis  were  aehieved  in  that  a  rather  eonvenient  and 
simplistie  model  was  created  to  expedite  the  process  to  create  an  initial  launch  trajectory 
for  the  first  stage.  The  optimal  trajectory  that  was  created  using  DIDO  was  deemed  to  be 
feasible  after  a  series  of  verification  and  validation  tests.  This  was  done  with  derived 
equations  from  Chapter  III  and  propagation  of  the  controls  and  plotting  the  results  against 
the  DIDO  solution.  Using  the  optimal  trajectory,  a  series  of  uncertainties  were  placed  on 
the  simulation  to  analyze  the  sensitivity  of  the  solution.  The  Monte  Carlo  simulations 
produced  small  deviations  in  the  endpoint  positions  and  had  little  effect  on  the  vehicle’s 
final  velocity.  From  here,  the  nominal  solution  could  be  sufficient  to  control  the  starting 
point  for  the  second  stage,  pending  additional  analysis.  There  is  substantial  room  for 
future  work  that  can  be  done  in  this  area. 

B.  HIGHER  FIDELITY  MODEL 

When  analyzing  the  dynamics  portion  of  this  model,  the  equations  used  were  a 
simplification  of  reality.  For  example  one  of  the  assumptions  in  the  model  was  that  the 
thrust  profile  is  constant  over  the  period  of  the  launch.  Depending  on  the  type  of  rocket, 
very  different  thrust  profiles  exist.  This  thesis  does  not  address  the  design  of  rocket 
motors  but  more  can  be  added  to  incorporate  thrust  profiles  for  various  rocket  motors 
used  in  industry.  More  emphasis  should  also  be  placed  on  the  aerodynamics  portion  to 
produce  more  accurate  approximation  of  aerodynamic  drag.  In  the  same  fashion,  models 
used  for  the  Earth’s  gravitational  force  did  not  account  for  the  oblateness  of  the  earth  or 
mass  distribution.  The  atmospheric  model  that  was  used  relies  on  a  sixth  order 
polynomial  created  as  a  curve  fit  base  on  observations  from  the  1976  standard 
atmospheric  model  that  only  goes  until  30,000  ft.  in  increments  of  1,000  ft.  While  the 
developed  model  seems  to  accurately  estimate  the  atmospheric  condition,  there  is  still 
room  for  increasing  the  fidelity  of  the  model.  While  most  of  the  time  higher  order  terms 
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are  neglected,  the  combination  of  multiple  “negligible”  terms  can  have  an  appreciable 
effect. 


C.  ADDITIONAL  STAGES 

Future  work  should  take  into  account  multiple  stages  to  create  a  more  realistic 
launch  to  orbit.  The  transition  from  single  stage  to  multiple  stages  would  be  rather 
seamless.  A  separate  function  block  would  need  to  be  written  such  that  the  starting 
position  for  the  subsequent  stages  would  be  the  ending  points  from  the  previous  stages.  In 
the  same  way  the  Monte  Carlo  simulation  could  be  run  to  assess  the  sensitivity  of  the 
control  for  each  stage. 

D,  CODE  ROBUSTNESS 

Other  launch  trajectory  generation  tools  should  be  compared  against  the  results  of 
this  thesis.  The  time  to  formulate  an  optimal  trajectory  would  be  the  most  important 
metric  to  determine  the  effectiveness  of  this  thesis  research.  Other  areas  might  include 
how  other  trajectory  generation  tools  perform  their  uncertainty  analysis. 
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